Ivermectin analysis
James Watson
3/11/2021
# Data were downloaded from http://diposit.ub.edu/dspace/handle/2445/173199
saint = read.csv("Ivermectin_study_data.csv", header=T)
# 1: ivermectin ; 2: placebo - we set placebo to be 1
saint$treat[saint$treat==2]=0
# take mean value of N and E PCR
saint = dplyr::filter(saint, !(is.na(ct_n) & is.na(ct_e)))
saint$ct = rowMeans(saint[, c('ct_n','ct_e')], na.rm = T)
ind_fit = saint$day<=21
iv_dat1 = data.frame(t = saint$day[ind_fit],
Trt = as.factor(saint$treat)[ind_fit],
id = saint$studyno[ind_fit],
y = saint$ct[ind_fit])run the model
CT_threshold = 32
data_stan = iv_dat1[iv_dat1$t<10, ]
ind_obs = data_stan$y<CT_threshold
ind_cens = !ind_obs
data_stan$id = as.numeric(as.factor(data_stan$id))
data_stan$Trt = as.numeric(as.factor(data_stan$Trt))-1
data_1 = list(Nobs = sum(ind_obs),
Ncens = sum(ind_cens),
J = length(unique(data_stan$id)),
id_obs = data_stan$id[ind_obs],
t_obs = (data_stan$t-1)[ind_obs],
Trt_obs = data_stan$Trt[ind_obs],
y_obs = data_stan$y[ind_obs],
id_cens = data_stan$id[ind_cens],
t_cens = (data_stan$t-1)[ind_cens],
Trt_cens = data_stan$Trt[ind_cens],
y_cens = data_stan$y[ind_cens],
CT_threshold=CT_threshold,
prior_Trt_sd=.1,
prior_CT_intercept = 20,
prior_sigma_sd = 1,
prior_mean_slope = 2,
prior_tau_intercept_sd = 2,
prior_tau_slope_sd = 1,
prior_tau_intercept_mean = 3,
prior_tau_slope_mean = .5)
interim_fit1 = sampling(object = stan_lin_model,
data = data_1,
chains=4, iter=10^5,thin=10)##
## SAMPLING FOR MODEL 'f8a746a2461316d8be14aecf5c16c64f' NOW (CHAIN 1).
## Chain 1:
## Chain 1: Gradient evaluation took 5e-05 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.5 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1:
## Chain 1:
## Chain 1: Iteration: 1 / 100000 [ 0%] (Warmup)
## Chain 1: Iteration: 10000 / 100000 [ 10%] (Warmup)
## Chain 1: Iteration: 20000 / 100000 [ 20%] (Warmup)
## Chain 1: Iteration: 30000 / 100000 [ 30%] (Warmup)
## Chain 1: Iteration: 40000 / 100000 [ 40%] (Warmup)
## Chain 1: Iteration: 50000 / 100000 [ 50%] (Warmup)
## Chain 1: Iteration: 50001 / 100000 [ 50%] (Sampling)
## Chain 1: Iteration: 60000 / 100000 [ 60%] (Sampling)
## Chain 1: Iteration: 70000 / 100000 [ 70%] (Sampling)
## Chain 1: Iteration: 80000 / 100000 [ 80%] (Sampling)
## Chain 1: Iteration: 90000 / 100000 [ 90%] (Sampling)
## Chain 1: Iteration: 100000 / 100000 [100%] (Sampling)
## Chain 1:
## Chain 1: Elapsed Time: 17.8136 seconds (Warm-up)
## Chain 1: 21.2163 seconds (Sampling)
## Chain 1: 39.0299 seconds (Total)
## Chain 1:
##
## SAMPLING FOR MODEL 'f8a746a2461316d8be14aecf5c16c64f' NOW (CHAIN 2).
## Chain 2:
## Chain 2: Gradient evaluation took 2.3e-05 seconds
## Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.23 seconds.
## Chain 2: Adjust your expectations accordingly!
## Chain 2:
## Chain 2:
## Chain 2: Iteration: 1 / 100000 [ 0%] (Warmup)
## Chain 2: Iteration: 10000 / 100000 [ 10%] (Warmup)
## Chain 2: Iteration: 20000 / 100000 [ 20%] (Warmup)
## Chain 2: Iteration: 30000 / 100000 [ 30%] (Warmup)
## Chain 2: Iteration: 40000 / 100000 [ 40%] (Warmup)
## Chain 2: Iteration: 50000 / 100000 [ 50%] (Warmup)
## Chain 2: Iteration: 50001 / 100000 [ 50%] (Sampling)
## Chain 2: Iteration: 60000 / 100000 [ 60%] (Sampling)
## Chain 2: Iteration: 70000 / 100000 [ 70%] (Sampling)
## Chain 2: Iteration: 80000 / 100000 [ 80%] (Sampling)
## Chain 2: Iteration: 90000 / 100000 [ 90%] (Sampling)
## Chain 2: Iteration: 100000 / 100000 [100%] (Sampling)
## Chain 2:
## Chain 2: Elapsed Time: 18.438 seconds (Warm-up)
## Chain 2: 21.152 seconds (Sampling)
## Chain 2: 39.59 seconds (Total)
## Chain 2:
##
## SAMPLING FOR MODEL 'f8a746a2461316d8be14aecf5c16c64f' NOW (CHAIN 3).
## Chain 3:
## Chain 3: Gradient evaluation took 2.3e-05 seconds
## Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.23 seconds.
## Chain 3: Adjust your expectations accordingly!
## Chain 3:
## Chain 3:
## Chain 3: Iteration: 1 / 100000 [ 0%] (Warmup)
## Chain 3: Iteration: 10000 / 100000 [ 10%] (Warmup)
## Chain 3: Iteration: 20000 / 100000 [ 20%] (Warmup)
## Chain 3: Iteration: 30000 / 100000 [ 30%] (Warmup)
## Chain 3: Iteration: 40000 / 100000 [ 40%] (Warmup)
## Chain 3: Iteration: 50000 / 100000 [ 50%] (Warmup)
## Chain 3: Iteration: 50001 / 100000 [ 50%] (Sampling)
## Chain 3: Iteration: 60000 / 100000 [ 60%] (Sampling)
## Chain 3: Iteration: 70000 / 100000 [ 70%] (Sampling)
## Chain 3: Iteration: 80000 / 100000 [ 80%] (Sampling)
## Chain 3: Iteration: 90000 / 100000 [ 90%] (Sampling)
## Chain 3: Iteration: 100000 / 100000 [100%] (Sampling)
## Chain 3:
## Chain 3: Elapsed Time: 17.9827 seconds (Warm-up)
## Chain 3: 21.4706 seconds (Sampling)
## Chain 3: 39.4533 seconds (Total)
## Chain 3:
##
## SAMPLING FOR MODEL 'f8a746a2461316d8be14aecf5c16c64f' NOW (CHAIN 4).
## Chain 4:
## Chain 4: Gradient evaluation took 2.4e-05 seconds
## Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.24 seconds.
## Chain 4: Adjust your expectations accordingly!
## Chain 4:
## Chain 4:
## Chain 4: Iteration: 1 / 100000 [ 0%] (Warmup)
## Chain 4: Iteration: 10000 / 100000 [ 10%] (Warmup)
## Chain 4: Iteration: 20000 / 100000 [ 20%] (Warmup)
## Chain 4: Iteration: 30000 / 100000 [ 30%] (Warmup)
## Chain 4: Iteration: 40000 / 100000 [ 40%] (Warmup)
## Chain 4: Iteration: 50000 / 100000 [ 50%] (Warmup)
## Chain 4: Iteration: 50001 / 100000 [ 50%] (Sampling)
## Chain 4: Iteration: 60000 / 100000 [ 60%] (Sampling)
## Chain 4: Iteration: 70000 / 100000 [ 70%] (Sampling)
## Chain 4: Iteration: 80000 / 100000 [ 80%] (Sampling)
## Chain 4: Iteration: 90000 / 100000 [ 90%] (Sampling)
## Chain 4: Iteration: 100000 / 100000 [100%] (Sampling)
## Chain 4:
## Chain 4: Elapsed Time: 17.9634 seconds (Warm-up)
## Chain 4: 22.2008 seconds (Sampling)
## Chain 4: 40.1642 seconds (Total)
## Chain 4:out=extract(interim_fit1)par(las=1, bty='n', family='serif',
mfrow=c(1,2), cex.lab=1.5, cex.axis=1.5)
xs = jitter(iv_dat1$t, amount = .3)-1
plot(xs, iv_dat1$y, pch=20, xlab = 'Time from enrollment',
ylab='Cycle threshold', xaxt='n', ylim=rev(range(iv_dat1$y)))
axis(1, at=c(0,3,6,13,20))
abline(h=32, lty=2, lwd=3, v=7, col='red')
for(id in unique(iv_dat1$id)){
lines(xs[iv_dat1$id==id], iv_dat1$y[iv_dat1$id==id],lty=2)
}
mtext(text = 'A)', side = 3, cex=1.5, adj = 0)
hist(100*out$beta_Trt, xlab='Treatment effect (%)',
main='', breaks = 50, yaxt='n',ylab = '', freq = F,
border = NA, col=adjustcolor('grey',.9))
lines(-300:300, dnorm(c(-300:300), mean = 0, sd = 100*data_1$prior_Trt_sd), col='blue',lwd=3)
abline(v=0, lwd=3, col='red')
mtext(text = 'B)', side = 3, cex=1.5, adj = 0)
writeLines(sprintf('Probability that treatment effect>0 is %s', mean(out$beta_Trt > 0)))## Probability that treatment effect>0 is 0.74765par(mfrow=c(5,5), mar=c(3,3,1,0),las=1, family='serif', cex.lab=1.5, cex.axis=1.5)
for(id in 1:data_1$J){
indo = data_1$id_obs==id
indc = data_1$id_cens==id
plot(data_1$t_obs[indo],
data_1$y_obs[indo],pch=20, xlab='', ylab='',cex=2,
ylim= rev(range(saint$ct)), xlim = c(0,6), yaxt='n', xaxt='n')
points(data_1$t_cens[indc],
data_1$y_cens[indc],pch=18,cex=2,col='blue')
polygon(c(c(data_1$t_obs[indo],data_1$t_cens[indc]),
rev(c(data_1$t_obs[indo],data_1$t_cens[indc]))),
c(c(apply(out$mu_obs,2,quantile,.1)[indo],
apply(out$mu_cens,2,quantile,.1)[indc]),
rev(c(apply(out$mu_obs,2,quantile,.9)[indo],
apply(out$mu_cens,2,quantile,.9)[indc]))),
border = NA, col = adjustcolor('grey',.3))
lines(c(data_1$t_obs[indo],data_1$t_cens[indc]),
c(apply(out$mu_obs,2,mean)[indo],
apply(out$mu_cens,2,mean)[indc]),lwd=3, col=adjustcolor('grey',.5))
points(data_1$t_obs[indo],data_1$y_obs[indo],pch=20,cex=2)
points(data_1$t_cens[indc],data_1$y_cens[indc],pch=18,cex=2,col='blue')
abline(h=CT_threshold, lty=2)
if(id%%5 == 1){ axis(2, at = seq(15,40,by=5))}
if(id>19) {axis(1, at = c(0,3,6)) }
}
Plot the posterior summaries
par(mfrow=c(2,3), cex.lab=1.5, cex.axis=1.5)
hist(out$beta_0, xlab='Population intercept',
main='', breaks = 50, yaxt='n',ylab = '', freq = F,
border = NA, col=adjustcolor('grey',.9))
lines(seq(0,30,by=.1),
dnorm(seq(0,30,by=.1),
mean = data_1$prior_CT_intercept, sd = 3),
col='blue',lwd=3)
hist(out$beta_t, xlab='Population slope',
main='', breaks = 50, yaxt='n',ylab = '', freq = F,
border = NA, col=adjustcolor('grey',.9))
lines(seq(0,10,by=.1),
dnorm(seq(0,10,by=.1),
mean = data_1$prior_mean_slope, sd = 1),
col='blue',lwd=3)
hist(out$beta_Trt, xlab='Treatment effect',
main='', breaks = 50, yaxt='n',ylab = '', freq = F,
border = NA, col=adjustcolor('grey',.9))
lines(seq(-1,1,by=.01),
dnorm(seq(-1,1,by=.01),
mean = 0, sd = data_1$prior_Trt_sd),
col='blue',lwd=3)
hist(out$tau_0, xlab='Standard deviation: random intercept',
main='', breaks = 50, yaxt='n',ylab = '', freq = F,
border = NA, col=adjustcolor('grey',.9))
lines(seq(0,10,by=.1),
dnorm(seq(0,10,by=.1),
mean = data_1$prior_tau_intercept_mean,
sd = data_1$prior_tau_intercept_sd),
col='blue',lwd=3)
hist(out$tau_t, xlab='Standard deviation: random slope',
main='', breaks = 50, yaxt='n',ylab = '', freq = F,
border = NA, col=adjustcolor('grey',.9))
lines(seq(-1,3,by=.01),
dnorm(seq(-1,3,by=.01),
mean = data_1$prior_tau_slope_mean,
sd = data_1$prior_tau_slope_sd),
col='blue',lwd=3)
hist(out$sigma, xlab='Standard deviation: residual error',
main='', breaks = 50, yaxt='n',ylab = '', freq = F,
border = NA, col=adjustcolor('grey',.9))
lines(seq(0,10,by=.1),
dnorm(seq(0,10,by=.1),
mean = 1, sd = data_1$prior_sigma_sd),
col='blue',lwd=3)